#coding:utf-8
from dataRead import readFeature, readLetter
import pandas as pd
import csv
import numpy as np

## m,n：所有数据的行数和列数，numcancer: cancer的行数；sequences:原始蛋白质序列
##返回result.行名为特征名，列名为[d,'ratio','flag']
def selectD(numcancer, gap, sequences):
    featurePath = r"graphs_model/colname.csv"#r"graphs_model/colname.csv"
    colnames = readFeature(featurePath)
    m = len(sequences)
    n = len(colnames)
    colnames.append('label')
    labels = np.concatenate( ( np.zeros((numcancer,1)), np.ones((m - numcancer,1)) ) )
    matrix = pd.DataFrame(np.append(np.zeros((m,n)), labels, axis = 1), columns =colnames, dtype = int)

    i=0
    for sequence in sequences:
        for j in range(0,len(sequence)-gap-1):
            feature = sequence[j]+ sequence[j+gap+1]
            matrix.ix[[i],[feature]]= matrix.loc[i, feature]+1
        i = i+1

    print(matrix)
    result =pd.DataFrame(np.zeros((n,3)), index = matrix.columns.copy().drop('label'), columns= ['d','rate','flag'])
    for col in matrix.columns.copy().drop('label'):
        D = -1
        maxRate = -1
        flag = False

        candidates = matrix[col]
        for candidateD in set(candidates):
            if(candidateD >= 1):
                count0 = len (  matrix[(matrix[col]>= candidateD) & (matrix.label==0) ] ) #统计cancer中指定特征个数>= candidateD的序列数
                count1  = len (  matrix[(matrix[col]>= candidateD) & (matrix.label==1) ]  )
                if(count1 == 0 or count0 == 0): #某类中有效序列数为0....
                    print( "count:" + str(count0) + " "+str(count1))
                    flag = True
                    continue
                rate = max( (count0*1.0/numcancer)*1.0 / (count1*1.0 / (m-numcancer)),(count1*1.0/(m-numcancer)) / (count0*1.0/numcancer)*1.0 )
                if(rate > maxRate):
                    maxRate = rate
                    D = candidateD

        result.ix[col] = [D,maxRate, flag]
    print(result)
    return result

import math
def selectFeature(dataFrame):
    m = len(dataFrame)
    mean = dataFrame.rate.mean()
    std = math.sqrt(dataFrame.rate.var()*(m-1)*1.0/m)  #无偏估计

    dataFrame['d_stand'] = dataFrame.rate.map(lambda x: (x-mean)*1.0/std)
    low = mean-2*std
    high = mean + 2*std
    features = dataFrame[(dataFrame["d_stand"] >= high) | (dataFrame["d_stand"] <= low) ]  #如何处理出现count0or count1 = 0的特征？？ | (dataFrame.flag == 1)
    print("出现d值使得有效的正类或者负类序列数为1的特征数：")
    print(len(dataFrame[dataFrame.flag==1]))

    return list(features.index)

'''
path = r"antioxidant/testMe.txt"
sequences = readLetter(path)[1]
m = len(sequences)
n = 441
selectD(m,n,4, gap=0, sequences= sequences)
'''
if(__name__ =="__main__"):

    pathPositave = r"E:\20170420蛋白质计算\virion1.txt"
    pathNegative  = r"E:\20170420蛋白质计算\non-virion1.txt"
    sequences = readLetter(pathPositave)[1]
    sequenceNeg = readLetter(pathNegative)[1]
    numAntioxidant = len(sequences)
    numAll = numAntioxidant+ len(sequenceNeg)

    sequences.extend(sequenceNeg)  #合并两个序列
    dataFrame = selectD(numAntioxidant, gap=0,sequences= sequences)
    featureSelect = selectFeature(dataFrame)
    print(featureSelect)
    print("2 sigma 外选取的特征数："+str(len(featureSelect)))

    '''
    path = r'test/sequence.txt'
    sequences = readLetter(path)[1]
    dataFrame = selectD( 4, gap=0,sequences= sequences)

    featureSelect = selectFeature(dataFrame)
    print("2 sigma 外选取的特征数："+str(len(featureSelect)))
    '''